Calculation of the optical response of C$o and Na 8 using time-dependent density 

functional theory and local orbitals 
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We report on a general method for the calculation of the frequency-dependent optical response of 
clusters based upon time-dependent density functional theory (TDDFT). The implementation is 
done using explicit propagation in the time domain and a self-consistent program that uses a lin- 
ear combination of atomic orbitals (LCAO). Our actual calculations employ the SIESTA program, 
which is designed to be fast and accurate for large clusters. We use the adiabatic local density ap- 
proximation to account for exchange and correlation effects. Results are presented for the imaginary 
part of the linear polarizability, Q ! a(aj), and the dipole strength function, S(u>), of Ceo and Nas, 
compared to previous calculations and to experiment. We also show how to calculate the integrated 
frequency-dependent second order non-linear polarizability for the case of a step function electric 
field, 7si ep (u;), and present results for Ceo- 



I. INTRODUCTION 

Although density functional theory (DFT)Bi is a very 
successful theory for the ground state properties, the ex- 
cited states calculated within the Kohn-Sham scheme of- 
ten are much less successful in describing the optical re- 
sponse and the excitation spectra. The solution to this 
problem, in principle, is the extension of DFT to the 
time-dependent systems. It is interesting to note that 
the first calculation^ using TDDFT preceded any for- 
mal development and it relied heavily on the analogy 
with the time-dependent Hartree-Fock method. The first 
steps towards the formulation of TDDFT were done by 
Deb and GoshtlO who fbeused on potentials periodic in 
time, and by BartolottiraH r»dio focused on adiabatic pro- 
cesses. Runge and Grossa established the foundations 
of TDDFT for a generic form of the timepcLependent po- 
tential. TDDFT was further developedSO to acquire a 
structure that is very similar to that of the conventional 
DFT. A very interesting feature of TDDFT, that does 
not appear in DFT, is the dependence of the density func- 
tionals on the initial state. For more information about 
TDDFT the reader is advised to read thep-authoritative 
reviews of Gross, Ullrich, and GossmannjLII and Gross, 
Dobson, and PetersilkaEj 

The polarizability describes the distortion of the charge 
cloud caused by the application of an external field. It is 
one of the most important response functions because it 
is directly related to electron-electron interactions, and 
correlations. In addition, it determines the response to 
charged particles, and optical properties. A quantity of 
particular interest is the dipole strength function, S(u>), 
which is directly related to the frequency-dependent lin- 
ear polarizability, a(ui), by 



a(u>) 



S{u')duj' 



(1) 
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The dipole strength function, S(lo), is proportional to 
the photoabsorption cross section, o~(u>), measured by 
most experiments and, therefore, allows direct compar- 
ison with experiment. In addition, the integration of S 
over energy gives the number of electrons, N ei {J- sum 
rule ) i.e. 



dES{E)=Y,fi 



(3) 



By taking the imaginary part of Eq. p) we obtain 



where /$ are the oscillator strengths. This sum rule is 
very important because it provides an internal consis- 
tency test for the calculations, indicating the complete- 
ness and adequacy of the basis set used for the compu- 
tation of the optical response. 

Optical probes are some of the most successful exper- 
imental tools that allow access to the properties of clus- 
ters. Consequently, there are many calculations of the 
optical response of small atomic aggregates. In particu- 
lar, there exist several theoretical studies of the examples 
chosen here, Ceo and Nas. This allows us to calibrate 
the accuracy of our method in comparison with other 
computational schemes. One of the first ab initio calcu- 
lations of the dipole response of atomic clusters -within 
TDDFT was performed by Yabana and Bertsch,E3 who 
studied large sodium and lithium clusters, and the Ceo 
molecule using a real time ajid space approach. For small 
Na clusters, Vasiliev et al&3 calculated the absorption 
cross section using the time-dependent density functional 
response theory (TD-DFRT) developed by Casida.lij 

The purpose of this work is to propose a method that 
will have significant advantages for the calculation of the 
polarizability of large clusters, reducing considerably the 
computation time while retaining the desired accuracy. 
This paper is organized as follows: In section O, we de- 
scribe the method of calculation. In section [n_C| we 
give details about the way we solve the time-dependent 
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Kohn-Sham equation and briefly summarize other meth- 
ods available. In section III , we present an overview of 
relevant calculations and the results of our calculation 
for Ceo and Nas- We compare our results with other cal- 
culations and experiments. 



IV, we describe 



In section 

the calculation and present the results for the imaginary 
part of the integrated frequency-dependent second order 
non-linear polarizability for the case of a step function 
electric field, S>7 step (o;), for Geo- In section [v|, we give 
the conclusions. 



II. METHOD OF CALCULATION 
A. Electronic structure calculations 

Our method involves the description of the elec- 
tronic states using linear combination of atomic orbitals 
(LCAO). Because the size of the LCAO basis is small, 
compared with other usual choices like plane waves or 
real space grids, the TDDFT calculations can be done 
efficiently using the techniques-described below. Our 
scheme is based on the SIESTAEljIia code, which is used 
to compute the initial wavefunctions and the Hamilto- 
nian matrix for each time step. SIESTA is a general- 
purpose DFT code which uses a local basis, and has been 
specially optimized to deal with large systems. As such, 
it represents an ideal tool for treating large clusters. Core 
electrons are replaced by norm-conserving pseudopoten- 
tialal3 in the fully nonlocal Kleinman-BylanderEj form, 
and the basis set is a general and flexible linear combi- 
nation of numerical atomic orbitals (NAOs), constr ucted , 
from the eigenstates of the atomic pseudopotentialsiiSEil 
The NAOs are confined, being strictly zero beyond a cer- 
tain radius. In addition, the electron wavefunctions and 
density are projected onto a real space grid in order to 
calculate the Hartree and exchange-correlation potentials 
and their matrix elements. 

The use of confined NAOs is very important for the 
efficiency of the SIESTA code. With them, by exploit- 
ing the explicit sparseness of the Hamiltonian and den- 
sity matrices, the computational cost for the construction 
and storage of the Hamiltonian and the electronic den- 
sity can be made to scale linearly with the number of 
atoms, in the limit of large systems. Therefore, a consid- 
erable effort has been devoted to obtain orbital bases that 
would meet the standards of precision of conventional 
first-principles calculations, while keeping their range as 
small as possible. A simple scheme for the generation 
of transferable bases that satisfy both requirements was 
presented in Refs. [l7] and |22|. These bases, which we 
utilize in this work, have been successfully applied to 
study the ground state properties of very different sys- 
tems, ranging from insulators-to metals, and from bulk to 
surfaces and nanostructures.113 It is not obvious however 
that these confined basis sets will be also adequate for 
the TDDFT calculation of the optical response. In this 



paper we show that, at least for the two systems consid- 
ered, the optical absorption can be accurately calculated 
using basis of NAOs with reasonable confinement radii, 
and a moderate number of orbitals per atom. Our results 
are in good agreement with other TDDFT calculations 
using computationally more demanding basis sets. 

Our approach is to carry out the calculations in the 
time domain, explicitly evolving the wavefunctions. We 
consider a bounded system in a finite electric field, i.e. 
the Hamiltonian includes a perturbation AH = — E • x. 
For the linear response calculations in this paper we have 
set the value of this field to 0.01 eV/A. The system is 
solved for the ground state usiog, standard time indepen- 
dent density functional theory.Eij Then we switch off the 
electric field at time t = 0, and for every subsequent 
time step we propagate the occupied Kohn-Sham eigen- 
states by solving the time-dependent Kohn-Sham equa- 
tion (h = 1) 

a* 

i^=H % (4) 
where H is the time-dependent Hamiltonian given by 



H 



1 -y 2 + V ext (r,t)+ [ t^Adr' + V xc [p](*,-t)- 



r - v 



(•5) 



The calculation of the exchange-correlation potential is 
done using the adiabatic local density approximation 
(ALDA) where V xc takes the form 



V xe \p](r,t) = 5pt[r) 



= V, 



LDA 



\pt](r). 



(6) 



E xF A [Pt] is the exch 



;e-correlation energy of the homo- 
geneous electron gasS It is important to notice that the 
Vxc in ALDA is local both in time and space. For every 
time step we solve Eq. (^), and from the new wavefunc- 
tions we construct the new density matrix 



(7) 



where cf (i) are the coefficients of the occupied wavefunc- 
tions which correspond to the basis orbitals cf)^ (r) . p^ v (t) 
has to be calculated and stored for overlapping orbitals 
only. The electron density is then obtained by 



p(r,t)=Xy"(t) ^(r)^(r), 



(8) 



and used for the calculation of the Hamiltonian in the 
new cycle. 



B. Calculation of the polarizabilities 

For every time step we calculate the dipole moment 
D(t) of the electrons in the cluster. This defines the 
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response to all orders and the frequency dependent re- 
sponse is found by the Fourier transform 



N-l 



D(w) = J dte luJt " St T>(t). 



(9) 



In our case we Fourier transform the dipole moment only 
for t>0. It is necessary to include a damping factor S in 
order to perform the Fourier transform. This damping 
factor gives the minimum width of the peaks of the imag- 
inary part of the response. Physically, it can be regarded 
as an approximate way to account for broadening. To lin- 
ear order the polarizability is given by D(w) = a(w)E(w), 
so that 



$fa(u>) 



i 

E 



(10) 



where the field is given by E(t) — E 9(—t). After Fourier 
transforming the dipole moment we obtain the elements 
of the frequency-dependent polarizability tensor aij(u>). 
We repeat the calculation with the electric field along 
different axis unless the symmetry is high enough that 
this is not needed. The average linear polarizability is 
given by 



I, 



< a(u) >= -Tr{atij(w)} 



(11) 



The choice of the coordinate system does not affect the 
average polarizability because of the rotational invariance 
of the trace. 



C. Solution of the time-dependent Kohn-Sham 
equation 

Efficient solution of the time-dependent Kohn-Sham 
equation (Eq. (0)) is of particular interest because to- 
gether with the calculation of the Hamiltonian, they are 
the most time consuming parts of the calculation. In this 
section we describe our approach of solving Eq. (m, as 
well as other existing methods used for the same purpose. 

In the LCAO formalism Eq. (Eh takes the form 



(12) 



where S is the overlap matrix between the orbitals and c 
is the column of the coefficients of the local orbitals. The 
overlap matrix is fixed for a given atomic configuration, 
hence we have to calculate and invert it only once. 
The formal solution of Eq. ( |l2|) is 

c(t) = U(t, 0)c(0) = Texp f-i J S^H^dA c(0), 

(13) 

where T is the time ordering operator. The most elemen- 
tary solution is obtained by breaking the total evolution 
operator into evolution operators of small time durations 



[/"(*, 0)~ JJ U((n + l)At,nAt), 



n=0 



where At = 



and 



U(t + At, t) = exp (-iS^HtyAt) 



(14) 



(15) 



T to t is the total time that we allow the system to evolve. 
The differences among propagation schemes arise from 
the way the exponential in Eq. ( |l5| ) is approximated. 
In our approach, we approximate the exponential in Eq. 
( |l5| ) with the Crank-Nicholson operator.EEl The coeffi- 
cients between the steps n + 1 and n are related by the 
equation 



„n+i _ 



1 - iS- l H{t n ) 



At 
2 



l + iS^Hitn)^ 



(16) 



This method is unitary, strictly preserving the orthonor- 
mality of the states for an arbitrary time evolution. For 
time independent Hamiltonians it is also explicitly time 
reversal invariant, and exactly conserves energy. In prac- 
tice, with a suitable choice of At, the energy is satisfacto- 
rily conserved even when the Hamiltonian changes with 
time. For example, in the calculations described below, 
the drift of the total energy at the end of the simulation 
(^130 fs in both cases) was only AE to t / E to t ~3x 10~ 7 
for C 60 and, - 8 x 1CT 6 for Na 8 , after N Ceo ~ 6100 and 
N^ 0s ~ 2800 time steps, respectively. The larger energy 
drift in the case of Nag is attributed to the use of larger 
time step. The method is stable when AtAE max << 1, 
where AE max is the range of the eigenstates of 
We can increase the stability of the solution if we in- 
clude more terms of the expansion in the numerator and 
denominator of the Crank-Nicholson operator, i.e. 



A/ 



Us- 



-l H At)2 
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At\3 



At 
2 



2 ' 



(17) 



By including more terms in the expansion it is possible ei- 
ther to increase the time step preserving the accuracy, or 
to increase the accuracy of the dynamics and the energy 
conservation for a given time step. The main advantage 
of using a bigger time step is the saving of time because 
we have to calculate the Hamiltonian fewer times. The 
energy resolution will not be affected since it depends on 
the total time that we allow the system to evolve. 

The method presented in this work has many sipajlar- 
ities with that described by Yabana and BertschjLia the 
main difference being the use of an LCAO basis set in 
the present case. However, this is a key difference be- 
cause the size of the matrices used in the calculations 
is considerably smaller compared to other basis choices. 
In addition, our method has other advantages associated 
with the real time formulation of TDDFT. Only occu- 
pied states are used in the calculation, in contrast to the 
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perturbative approachll3e3 where there is a sum over the 
excited states of the system. The implementation is rela- 
tively simple, since we use essentially the same operations 
as already used to find the ground state properties. It is 
also advantageous that nonlinear effects can be included 
in a straightforward way. One disadvantage of the real 
time approach is the calculation of the Hamiltonian for 
every time step. Although this is not an attractive fea- 
ture there is no other way to calculate the time evolution 
of the system. 

There are many other ways to approach the solution 
of the equations. For completeness, we discuss in the 
Appendix several methods that could be of potential rel- 
evance for our calculations. One of the main goals of 
these methods is to enable longer time steps. This would 
be a great advantage in our work; however, the fact that 
the self-consistent Hamiltonian changes as a function of 
time, limits their use. 
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FIG. 1. Dipole strength function of Nag vs. energy. 
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III. DISCUSSION OF RESULTS 
A. Small metal clusters: Nag 

The first calculation we performed is the optical re- 
sponse of Nag. The main purpose of this calculation was 
to investigate the accuracy of our method in the case of 
a small cluster, where the effects related to the confine- 
ment of the orbitals should be more noticeable and where 
the size of our basis is much smalles than in previous cal- 
culations using real space grids .til Na§ is the smallest 
closed shell Na cluster that its optical response exhibits 
the presence of a_p]asmon which is experimentally ob- 
served at 2.53 eV£3p The width of the plasmon is due 
to Landau damping. E2 

Previous work can be grouped into two types: earlier 
work on jellium spheresE3'E2l~El that reproduces the qual- 
itative features but not the .quantitative energies of the 
peaks, and more recent workliJS that takes into account 
the detailed atomic structure and— }S-4n general in very 
good agreement with experiment. OS IPnf^e nrs t cate- 
gory are the calculations of Selby et aLj2Zll23 who calcu- 
lated the photoabsorption cross section using the modi- 
fied Mie theory, which is a classical theory. The plasmon 
was found to be at ~ 2.76 eV. By using the self-consistent 
jellium model in the time dependent local density .approx- 
imation (TDLDA), first introduced by EkardtjHj Yan- 
nouleas et alM calculated the photoabsorption cross sec- 
tion of Nas and predicted the plasmon position at 2.82 
eV. 

Bounacic-Koutecky et al^ calculated the absorp- 
tion spectrum of Nas using the configuration-interaction 
method (CI). Because the all-electron calculation is com- 
putationally very demanding, they obtained the excited 
states by a non-empirical effective core potential cor- 
rected for the core- valence correlation using a core po- 
larization potential. The position of the plasmon was 
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FIG. 2. Imaginary part of the linear polarizability of Nag 
vs. energy. 

predicted at ~ 2.55 eV. Vasiliev et aL0 calculated the 
photoabsorption cross section using TD-DFRT.E3 Their 
calculations made use of norm-conserving pseudopoten- 
tials and a real-space mesh as a basis set. The position 
of the plasmon agreed, with the photoabsprption experi- 
ments of Selby et alH and Wang et alJB within 0.1-0.2 
eV. 

In our calculation we let the system evolve for the 
total time of T=31.42 eV" 1 . The energy resolution, 
determined by Alu — ^, is, in consequence, equal 
to 0.1 eV. The time step is 11.025 xl0~ 3 eV _1 , and 
the damping factor used in the Fourier transform is 
0.095 eV. Troullier-Martins pseudopotentialsEj including 
non-linear partial core corrections^ for the exchange- 
correlation interaction between valence_and core elec- 
trons, and an auxiliary real-space gridEj equivalent to 
a plane-wave cutoff of 70 Ry are also used in this 
calculation. The basis set includes 13 NAOs per 
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atom: two radial shapes to represent the 3s states 
plus a polarizationO p shell with confinement radii 
r s —rp° L — 12.2 a.u., and two additional 3p and 3d shells 
with radii r p —rd— 10.0 a.u.. 

Fig. 1 and Fig. 2 present, respectively, our results for 
the dipole strength function and the imaginary part of 
the linear polarizabihty of Nas for energies up to 4 eV. 
The shape of these curves is in excellent agreement with 
both, the calculations_pf Vasiliev et alt3 and, the experi- 
ments of Wang et al.L3 However, the results appear to be 
shifted to higher energies. In fact, the maximum of the 
plasmon peak is obtained at 2.86 eV, which is 0.33 eV 
higher than the experimentally observed value. This shift 
to higher energies seems to be related to the extension 
of the LCAO basis: using more confined orbitals we get 
a larger shift. The integrated dipole strength is equal 
to 6.97 out of 8, thus fulfilling 87.13 % of the sum rule. 
The partial fulfillment of the sum rule signifies the incom- 
pleteness of our basis set. The static linear polarizabihty 
a(0) can be obtained from standard (static) calculations 
of the induced dipole as a function of the applied field. 
Using this approach we obtain a value of 13.2 A 3 /atom. 
An alternative way to calculate a(0) is provided by the 
formula 

^^r^^r^, as) 

m Jo w T Jo w 

from which we obtain a value of 12.5 A 3 /atom. (This 
result can also be derived from the fact that for the step 
perturbation D(t = 0) = a(0)E.) The discrepancy be- 
tween both estimations is probably related with the lack 
of energy resolution of the calculated a(uj) to perform the 
integral in Eq. ( |is| ) with the required accuracy. Both 
results are in reasonable agreement with the values of 
14.6 l—Zatom and 14.7 A 3 /atom, computed by Vasiliev 
et alMM using the TDLDA and finite field methods, re- 
spectively. 

B. Large molecules: Ceo 

The best known Buckyball Ceo is a very interesting 
system with strong electron-electron interactions due to 
the confinement. There are quite a few calculations con- 
cerning the optical properties of Cgo and in particular 
its optical response. The main feature of the optical re- 
sponse of Ceo is the presence of two collective excitations 
(plasmons). The low energy plasmon can be associated 
with the 7r electrons while the high energy plasmon with 
both the a and ar electrons, in analogy with the plasmons 
in graphite.BEj The plasmeas have been observed by a 
plethora of experiments.EZTEil .— . .— . 

The earliest theoretical workc3c3 on Ceo involved sim- 
plifying approximations for the electron states (tight- 
binding or spherical averaging) and for the electron in- 
teraction (neglect or RPA-like treatmentsl—We will com- 
pare our results with those of Westin et a/O and Yabana 



and Bertschliia, who used large basis sets and realistic 
carbon potentials. Westin et al. used single particle 
wavefuctions, determined from a local density approx- 
imation (LDA) calculation, to evaluate the dipole ma- 
trix elements which combined with a sum over states ap- 
proach yielded the unscreened frequency-dependent lin- 
ear polarizabihty. Screening was included in a RPA-like 
fashion by introducing an effective screening parameter. 
The polarizabihty calculated in the static limit was used 
to evaluate this parameter for the calculation of the dy- 
namic response. The optical response and the sum rule 
for the low energy part were in reasonable agreement with, 
the experiment of Leach et al.EJ Yabana and BertschE 3 ] 
used TDLDA, evolving the system in real space and time, 
to calculate the dipole strength function of Ceo- Their 
calculation also gives reasonable-agreement with the ex- 
perimental data of Leach et a/.Jiil for the sum rule of the 
low energy part although it misses many details of the 
structure. This calculation is very similar in quality to 
ours. 

The total simulation time in our calculation of the po- 
larizabihty and dipole strength of the Cgo molecule is 
again 31.416 eV -1 , and the corresponding maximum en- 
ergy resolution of 0.1 eV. The time step however, which 
is set equal to 5.145 xl0~ 3 eV -1 , is smaller than the one 
used for Na§. This is because of the higher frequency 
range of the response of Cgo- The damping factor used 
in the Fourier transform is equal to-0.34 eV in this case. 
Troullier-Martins pseudopotentialstJ, a double-^ polar- 
ized basis set, and a real-space grid cutofO of 70 Ry 
were used in this calculation. There are 13 NAOs per 
C atom: two different radial shapes for the description 
of the 2s states, another two for the 2p, plus an addi- 
tional shell of d orbitals. The radii of confinement used 
are r s =5.12 a.u. and r p =r^ oi - =6.25 a.u. (corresponding 
to an energy shifvQ of 50 meV). For Cgo, the calculated 
spectra show small dependence in these radii, at least as 
far as they are not selected to be very stringent. In Fig. 3, 
the dipole moment is shown as a function of the time step 
number. The dipole strength function obtained from the 
time evolution of the dipole moment is shown in Fig. 4 
for energies up to 60 eV. Its main features are the low 
energy transitions that come from the it electrons and 
the a and ir electron transitions in the region of 14-27 
eV. In the low energy part of dipole strength function we 
have peaks at 3.46, 4.35, 5.36, and 5.84 eV, which agree 
very well with-the ones obtained by the calculations of 
Westin et al.L3 By integrating the dipole strength func- 
tion over energy we get the sum rule strength. The total 
sum rule strength is 223.78 out of 240. Therefore, we 
satisfy the sum rule up to 93.24 %. This reflects the in- 
completeness of our basis set, which fails to reproduce 
some of the excitations in the high energy part of the 
spectrum. The a plasmon is broadened, but this is a 
common feature of all the TDDFT calculations done for 
C60-H3 In Fig. 5, the imaginary part of the polarizabihty is 
given as function of energy. By using Eq. (|l8|), the static 
linear polarizabihty a(0) is found to be 91.1 A 3 , while 
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FIG. 3. Dipole moment of Ceo vs. number of time steps. 
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FIG. 4. Dipole strength function of Ceo vs. energy. 
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FIG. 5. Imaginary part of the linear polarizability of Ceo 
vs. energy. 



our finite field calculations produce a value of 87.3 A 3 . 
Results for a(0), from very accurate finite- field calcula- 
tions using fifteen values of the field, are given in section 
[Tv| . Both values are higher than the lower limit estima- 
tion of 62.5 A 3 from quantum- mechanical calculations ,cZl 
and in good agreement with the value of 85 A 3 obtained 
by Yabana and Bertsch.EJ They also agree well with the. 
experimental values of 79.3 A 3 rfetfl UV absorption,E3 
and 85.2 A 3 from EELS spectra.EEJ 



IV. NON-LINEAR POLARIZABILITIES 

Because of the non-perturbative nature of our method 
we are able, for large values of the applied field, to obtain 
non-linear polarizabilities. In this section, we present 
the calculation of the imaginary part of the integrated 
frequency-dependent second order non-linear polarizabil- 
ity, Cry ste p(w), which is related to the response to a step 
function, for Ceo- Because Ceo is centrosymmetric the 
first order non- linear polarizability, [3{uj), and all other 
polarizabilities involving an even number of fields, van- 
ish by symmetry. 

The advantage of the explicit time method is that ex- 
actly the same methods can be used to calculate the non- 
linear response of the system. The disadvantage is that 
(unlike the linear case where each Fourier component is 
independent) the non-linear response depends upon the 
detailed spectrum of the applied field. Here we derive 
the non-linear response of an electric field coupled to a 
Ceo for the case where the field is the step function used 
before. A different calculation would have to be done 
to find the non-linear response to a field with a different 
time dependence. 

First we give the relation of our calculation to the gen- 
eral definition of second order non-linear response, as a 
function of time, which iscil 



D^{t) = f dti f dt 2 f dt 3 
1 {t-t l MM)E(t 1 )E(t 2 )E{t 3 ). 



(19) 



For the case of a step function perturbation i.e. E(t) = 
E 0(—t), it takes the form 



L> (3) (i) = iE 3 lim 



(20) 



(kiJxdu)2(kiJ3 

s t -o+ J (2tt) 3 
e -^ 1 +^ 2 +^ 3 )* 7 (_^ 1 _ U2 _ u ljU j 2t 0,3) 

(wi - i5i)(u> 2 - - iS 3 ) 



We Fourier transform Eq. ( |20| ) and obtain the second 
order non-linear response as a function of frequency 



iE 3 lim 



s z ^a+ J (2tt) 2 
7(-w; ui -u> 2 - uj 3: uj 2 , W3) 
(u — uj 2 ~ 1V3 - i5i)(uj 2 ~ iS 2 ){uj3 — 163) 



(21) 



G 



The quantity we calculate is the real part of the second 
order non-linear response, sRD^iuj), from which we can 
extract the imaginary part of the integrated second order 
non-linear polarizability, ^sj step (uj). Explicit details are 
given below and in analogy to Eq. ( |i"o| ) 37 step (w) is given 
by 



1 .5e-35 



-u) lim 



57(-cj; uj - u 2 - u>3, lo 2 , UJ3) 
(bj - UJ2 — 0J3 — iSi)(u>2 — iS 2 )(uJ3 — 163) 



du}2<lu}3 



(22)! 



Just as in Eq. ( fL8| ) for the linear term, 3 ! 7 s t ep (a;) can be 
related to the static second order non-linear polarizability 
by the expression 
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%lste P (to) =7(0; 0,0,0). 



(23) 



Eq. (|23|) can be trivially derived by realizing that 
D^ 3 \t = 0) = 7(0; 0, 0, O)^ 3 when a step function pertur- 
bation is applied. Alternatively, we can derive Eq. j2^ ) 
directly from Eq. ( pl| ) by applying the Kramers-Kroning 
relations for j(—u)i — 0J2 — 0J3; u>i, UJ2, W3). In fact, with 
the help of the Kramers-Kroning relations we can derive 
another interesting equality for the first moment of our 
integrated response, 



1 



dio ^sj s tep(u) = du 3ty(— oj;oj,0,0). 



(24) 



For the calculation of 7 s tep(w) we calculate the re- 
sponse of the system under two different step function 
perturbations. In the first calculation the field used is 
equal to Ei = 0.10 V/A, and we assume that the re- 
sponse Di(u>) is linear with respect to the field. This 
assumption is true since for to = the non-linear terms 
contribute to the response only 4.82xl0 _3 %. The con- 
tribution is of the same order of magnitude for u> ^ 0. 
In the second calculation the field is equal to E% = 1.00 
V/A, and we assume that the response D 2 (uj) consists of 
the linear response and the second order non-linear re- 
sponse D 2 3 \lu). The values of the field are in the same 
range as those used by Westin et alsB for the determina- 
tion of the static second order non-linear polarizability. 
Using Eq. (|i~0|), we have that 



and 



Ei 

D\(uS) — a(uj)Ei(uj) = a(ijj)- — 



D 2 {lo) - a(w)E 2 (u) = \u) 



From Eq. <M) it follows that 



D 2 3 \co) 



E. 



2 I 



(25) 



(26) 



(27) 



and from Eqs. (pa), (Eq), and (p7|) we obtain 
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FIG. 6. Imaginary part of the integrated second order 
non-linear polarizability , S?7 s tep(w), of C60 vs. energy. 
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(28) 



Our calculation for 7 s t ep (w) is quite straightforward in 
contrast to the perturbative method where it becomes 
computationally very demanding. 

In Fig. 6, we present the results, up to 60 eV, for 
97ste P (w), where Ss7 step (a;) is given by Eq. (H|). As ex- 
pected, 97 s tep(w) has both positive and negative values. 
The reason why ^s^ s t ep (uj) does not vanish below some 
finite frequency (as does the linear response) is because 
the second-order non-linear term represents many pro- 
cesses of both absorption and emission of photons and 
the Ceo molecule can couple to a continuum of modes 
that extends to zero frequency. This can also be seen in 
the integral expression Eq. p2]). 

Similarly to the case of the linear polarizability, we 
can obtain an estimation of the magnitude of 7(0; 0, 0, 0) 
from static self-consistent calculations performed with fi- 
nite fields. This value can be contrasted to similar calcu- 
lations in the literature, providing an estimation of the 
uncertainty of our calculations of the non-linear terms, 
and an internal consistency test for our calculation of the 
integrated frequency-dependent second order non-linear 
polarizability. We have followed here a procedure similar 
to that used in Ref. |5^, performing LDA calculations of 
the total energy and electric dipole of the Cgrj molecule 
for fifteen different values of an external static electric 
field E ranging from 0.003 V/A to 2.0 V/A. The results 
of the total energy were then fitted using the expression 
Wtot =Wo-^aE 2 -|;7E 4 , where a is the linear polarizabil- 
ity and 7 the second order non-linear polarizability. The 
values obtained for a and 7 are, respectively, 85.3 A 3 and 
3.54xl0 -36 esu. The results for the dipole moment were 
fitted to the expression D=aE+7E 3 . The corresponding 
values obtained are 85.3 A 3 and 3.71 xlO -36 esu. These 
values of 7(0; 0, 0, 0) are in excellent agreement with the 
value of 3.72xl0~ 36 esu obtained using Eq. (23), thus 
confirming the validity of our calculation for 7 s t ep (w). 
All of the values obtained by our calculations are smaller 
than the upper bound of 3.7xl0~ 35 esu proposed by the 
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experiments of Geng and Wright.E3 Also, our results for 
7(0; 0, 0, 0) are in reasonably good agreement withi-ether 
LDA calculations in the literature. Quong et al£3 re- 
ported values of 82.7 A 3 and 7.0xl0~ 36 esu, for a and 7, 
respectively, using an all-electron method with a Gaua= 
sian expansion as a basis set. Van Gisbergen et al£3 
reported very similar values, 82.5 A 3 and 7.3 xlO -36 esu, 
using a computational scheme based on a frozen-core ap- 
proximation, and a basis set of Slater functions. It is 
interesting to note that these results are much smaller 
than those obtained using simplified tight-binding mod- 
els within an independent electron picture, where the 
effects of screening are neglected. In such calculations 
the values of yiQ'A, 0, 0) obtained are of the order of 



200x10" 



esu 



V. CONCLUSION 

We presented a method for the calculation of the op- 
tical response of atoms and clusters. The main features 
of the method are the description of the wavefunctions 
in terms of an efficient local orbital (LCAO) basis and 
the explicit evolution of the system in time. This ap- 
proach is designed for large clusters and in fact it gives 
excellent results for Ceo- It is also shown to work remark- 
ably well even for systems for which the LCAO basis is 
very small, such as Nas- Our approach has the desirable 
features that only occupied states are needed and that 
all the most computationally intensive operations are es- 
sentially the same as those used to calculate the ground 
state properties. In addition, non-linear effects can be 
included in a straightforward way, and we have shown 
how to calculate the second order non-linear response for 

C60- 
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APPENDIX: ALTERNATIVE APPROACHES FOR 
THE SOLUTION OF THE TIME-DEPENDENT 
KOHN-SHAM EQUATION 

Other ways to approximate the exponential in Eq. ( |l4| ) 
include the expansion p£ the exponential in a series of 
Chebyshev polynomials!^ 



exp (-iS~ l HAt) ~ exp (-i(AE/2 + E min )At) x (Al) 

71 = ^ ' 

where <p n are the Chebyshev polynomials, and the ex- 
pansion coefficients a n (x) can be shown to be analo- 
gous to Bessel functions of the first kind of order n. 
Hnorm is a normalized Hamiltonian 2(S~ 1 H — E av )/ AE, 
where E av — (E max -t- E nl ^ n ) AE — E max E rn ^ n . E max 
and E m i n are, respectively, the maximum and minimum 
eigenvalues of S^ 1 !!. The Chebyshev polynomials are 
chosen because their error decreases exponentially when 
N is large enough, duetto the uniform character of the 
Chebyshev expansion.Ej Time reversal is built into the 
expansion coefficients, but the method does not effec- 
tively conserve norm or energy. While it wocks remark- 
ably well for time- independent Hamiltonians,EJ for time- 
dependent Hamiltonians the method becomes inefficient 
as the Chebyshev polynomials of the Hamiltonian have 
to be recalculated for each time step. 

Another approximation for the exponential in Eq. ( |T4| ) 
is performed by using the split-operator method intro- 
duced by Feit et al!Ei According to this method, the ex- 
ponential which contains the Hamiltonian operator can 
be split as 



exp[-i(T + V)At] ~ 

exp(— i— TAt) exp(-iVAt) exp( 



(A2) 



-i\TAt). 



This method was later generalized by SuzukO for an ar- 
bitrary number of oper ators, providing higher order ex- 
pansions (formula (A2) corresponds to second order in 
At), and a rigorous extension of the method to time- 
dependent Hamiltonians.L^I The split operator method 
takes advantage of the fact that it is very convenient to 
treat operators in their diagonal representations. For ex- 
ample, it is trivial to apply the kinetic energy operator 
to a wave-function in Fourier space, while the effect of 
a local potential is more easily calculated in real space. 
This method is, in principle, unconditionally stable and 
norm conserving. On the other hand, it does not con- 
serve energy and it can only be used to Hamiltonian op- 
erators which can be split into two non-commuting parts 
with a simple transformation between them. Therefore, 
the method is very well suited for plane wave or real 
space methods, where efficient fast Fourier transform al- 
gorithms provide an exact transformation between finite 
plane- wave expansions and real space grids. In other 
words, they span exactly the same subspace of functions, 
and it is possible to switch between representations where 
the kinetic energy and the potential are diagonal within 
a given subspace. For an LCAO basis this is not possi- 
ble. If we take an arbitrary wavefunction expanded in an 
orbital basis tp{r) = ^ v (py^), the result of applying 
the operator exp(-iVt) using a grid of points in real- 
space /(rj) = exp(— iV(ri)t)ip(ri) will be, in general, not 
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representable using the same local basis, i.e. some of 
the resulting function has been spilled from the subspace 
spanned by the local basis. 

An alternative way to solve Eq. (|lj) is by using the 
second-order differencing (SOD) method introduced by 
Askar and CakmakS In SOD the symmetric relation is 
used 



z{t + At) - c(t - At) = (e" 



-iS~ 1 H(t)At _ iS^ 1 H{i)At 



)c(t) 
(A3) 



and by expanding the exponential in Taylor series, the 
second-order propagation scheme is obtained 



c(t + At) ~ c(t - At) - 2iAtS~ 1 H(t)c{t). 



(A4) 



By construction SOD obeys the time-reversal symmetry. 
SOD is stable, and norm and energy are approximately 
conserved, only if the Hamiltonian is Hermitian. In spite 
of its simplicity, an important drawback of this method 
for its application to large systems is that it is necessary 
to store the wavefunctions at two different time steps, 
which may require a large amount of memory. 

Finally, another method, for solving Eq. (|l^) is the 
short iterative Lanczos.EJE-3 It is very convenient for cal- 
culations that involve very big Hamiltonians, especially 
when they are time independent. The Hamiltonian is 
projected to a subspace of smaller dimensionality and 
takes a tridiagonal form that makes easy to perform cal- 
culations. The Lanczos recurrence creates a set of orthog- 
onal polynomials which constitute a finite polynomial ap- 
proximation of the operator. An interesting feature of the 
method is its dependence on both the operator and the 
initial vector. The Lanczos recurrence relation is 



(S 1 H) qj = ^-iqj-x + aj-qj + ^qj+i- 



(A5) 



The coefficients are aj — (qj, (S _1 if) qj) and /3,-_i = 
(qj— l, (S~ 1 H) qj), where (a, b) = Sb is the usual com- 
plex inner product, and (qi,qj) = <5y. Matrices H and 
S have Nb x Nb dimension and the vectors qj have Nb 
components, Nb the total number of basis functions. For 
very large Nb, where the method is particularly useful, 
the inverse S^ 1 is not directly calculated but an itera- 
tive method is used instead.O The recurrence relation is 
initiated by setting 



and 



qo = c(0) 



(S 1 H) q = a q + /?oqi- 



(A6) 



(A7) 



After P iterations, the projected operator (5 _1 iJ)p = 
(qj, (S 1-1 H) qi) is a P x P matrix with a tridiagonal 
form that can be very easily diagonalized. The solution 
of Eq. dl2|) becomes 



where Z is the NbX P transformation matrix that diago- 
nalizcs (S~ 1 H)p and Dp is the diagonalized matrix. For 
an arbitrary initial condition, the accuracy of the time 
evolution achieved with Eq. (AS) is equivalent to a P 
order expansion of the evolution operator 



p-i 



c(At) = J2 ~ir( s ' lH ) k c (°)> 



k=0 



k\ 



(A9) 



and, consequently, the energy is only approximately con- 
served. However, the propagator in Eq. (A8) is unitary, 
and the normalization condition is strictly preserved. For 



time independent Hamiltonians, Eq. (A8) can be used to 



evolve the wavefunction for a time interval r that depends 
on P. For times larger than r, the time evolution pre- 
dicted by a P order expansion becomes inaccurate, and 
it is necessary to recalculate propagator using c(r) as the 
starting point for the recurrence procedure. The method 
is, therefore, very well suited to follow the evolution of 
wavefunctions described by a large basis set during short 
periods of time (< r). 

For the purpose of the calculation of response func- 
tions, we believe that the method adopted in this paper, 
based on the Crank-Nicholson operator, is superior to the 
short iterative Lanczos, at least for moderate basis sizes. 
There are several reasons for this: i) We need to evolve 
the wavefunctions for long times, ii) The Hamiltonian is 
time dependent. This implies that the recurrence cycle 
for the calculation of the approximated time evolution 
operator has to be repeated for each time step. Hi) In 
our case, we have to evolve all the occupied states. The 
propagator in Eq.(A8), however, depends on the initial 
state. Therefore, it is necessary to develop some general- 
ization of the scheme presented. For example, it might be 
also possible to construct an approximate time evolution 
operator starting from some weighted linear combination 
of the occupied states, and projecting it into a subspace 
of dimension P > N oec , where N occ is the number of 
occupied wavefunctions. 



c(At) = Z T e 



MO), 



(A8) 
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